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Abstract. We study a Lagrangian numerical scheme for solution of a nonlinear drift diffusion 
equation on an interval. The discretization is based on the equation's gradient flow structure 
with respect to the Wasserstein distance. The scheme inherits various properties from the 
continuous flow, like entropy monotonicity, mass preservation, metric contraction and minimum/ 
maximum principles. As the main result, we give a proof of convergence in the limit of vanishing 
mesh size under a CFL-type condition. We also present results from numerical experiments. 



1. Introduction 

In this paper, we propose and analyze a very particular spatio-temporal discretization of the 
following nonlinear initial-boundary value problem on an interval / — [a, b]: 

dtU^P{u),, + {V,{x)u),, u.,{t;a)=u,it;b)^0, u{0;x) ^ u°{x) > 0. (1) 

Specifically, we are interested in approximating non-negative weak solutions w : [0, T] x / — > M>q 
to ([T]) on arbitrary time horizonts T > 0. Our assumptions are that 

• the nonlinearity P : M>o M is continuous, is C^-smooth on M>o, and satisfies 
P(0)=0, P'(r) > 0, limP'(r)<oo, lim P'(r) = +oo, si-^P(l/s) is concave, (2) 

the prototypical example being the porous medium term P(r) ~ r'" with some to > 1; 

• the potential y : / — > M is C^-smooth with 

V^{a)^V,.{b)^0. (3) 

Under the given regularity assumptions, there are numerous possibilities to design efficient numer- 
ical schemes for solution of ([T]), e.g., using finite differences. The particular discretization under 
consideration here is special insofar as it is based on the representation of ([T]) as a gradient flow, 
and it inherits certain qualitative features of that variational structure; see below. 

1.1. Gradient flow structure. We summarize some basic facts about the variational formulation 
of 0. The divergence form in combination with the no-flux boundary conditions and ([3| implies 
the conservation of mass 



J u{t; x) dx = Af := J u°{x) dx > for aU t > 0, 



(4) 



and we shall consider M as some fixed quantity from now on. The space D^^(/) C L^{I) of bounded 
and strictly positive densities u with total mass M can be endowed with the L^-Wasserstcin metric 
W2; the definition and elementary properties of this metric are reviewed in Section[2]below. Next, 
introduce the energy functional E for u G D*^(/) by 



E(u) = j (t){u{x))dx + J u{x)V{x)dx, (5) 

where the internal energy potential : M>o — > M is an arbitrary second anti-derivative of r h- > 
P'(r)/r; see 

The link between the energy E and equation ([T]) — which has been rigorously established 
by Otto [E] — is that solutions to ^ form a gradient flow in the energy landscape of E with 
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respect to the metric W2. Further, it has been observed by McCann [IT] that the functional E is 
(— A)-convex along geodesies in W2, with 

A = max(-Kx(2;)) >0. (6) 

Consequently, the gradient flow is (— A)-contractive. Some implications are: 

(1) The energy E{u{t)) is monotonically decreasing in t. 

(2) Two solutions w, v diverge at most at an exponential rate of A in the Wasserstein distance, 
i.e., 

W2(u(i),i;(t)) < W2(u°,w°)e^* forant>0. (7) 

(3) There is a unique non-negative and mass preserving global solution for measure valued 
initial conditions, i.e., the density w*^ in ([!]) can be replaced by an arbitrary non-negative 
measure on / with mass M. 

Below, we discuss in which sense these properties are inherited by our discretization. 

1.2. Discretization. Semi-discretization in time of gradient flow equations like ^ has become 
a key tool in existence proofs and for the rigorous derivation of a priori estimates. The celebrated 
minimizing movement scheme [T] (also referred to as JKO |13l or simply implicit Euler scheme) 
works in the situation at hand as follows: given a time step r > 0, one defines inductively — 
starting from m° = — approximations it" of u{nT) as minimizcrs in D*^(/) of "penalized 
energy functionals" £,-(•, given by 

E,K<-i) = ^W2K<-i)2+E(w). (8) 

Thanks to the (— A)-convexity of E, it follows from the theory developed in [1] that the functions 
Ur : [0,00) — >■ D*^(/) obtained by piecewise constant interpolation in time converge for r | to 
the unique weak solution u : [0, 00) — > 'D'^{I) of Q. 

To obtain a full (spatio-temporal) discretization, we perform the minimization of E^- not over 
the entire set D*^(/), but over a submanifold D^-'^(/) of finite dimension [K — 1) G N. The 
discretization parameter ^ = (Coi ^ii • • • j Cif) is an increasing sequence of numbers G [0,M], 
with = ^0 < ?i < ^2 < • ■ ■ < C/f = The corresponding submanifold D^^(/) consists of 
piecewise constant density functions u G D^^(/) of the form 

K 
k=l 

where the end points Xk of the intervals are variable subject to the constraint 

a — xq < xi < ■ ■ ■ < xji ~ 6, 

and the positive weights Uk are given in terms of the x^ by 

{xk - Xk-i)uk = 4 := Cfc - Cfe-i > 0. 

One may think of Si to Sk as lumps of mass, each of which is uniformly distributed on its respective 
interval {xk-i,Xk], see Figure [l] We refer to this discretization as Lagrangian scheme. 

Given a discretization A = consisting of a time step r > and a spatial mesh ^, and 

an initial condition vP^ G D|^(/), define a discrete solution ua = ('^A' ^A' ■ ■ ■) inductively by 

= argmin 'Eriu, u'^^) for n = 1, 2, . . . (9) 

«GDf(/) 

The seemingly involved definition of the recursion Q leads to a simple and practical numerical 
scheme, whose complexity is comparable to that of a standard discretization of ([T]) by finite 
differences. In addition, the Lagrangian nature of the scheme admits a geometric interpretation 
of the solution in terms of transportation of mass elements on / along the characteristics Xk(t). 
Concerning structure preservation, we summarize some noteworthy features of this approach: 

• The energy ^{u^) is monotone in n, and all are positive and have the same mass M. 
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• Any two discrete solutions and va on the same grid A satisfy the contraction estimate 

^2{ul,vl) < (l-2Ar)-"/2W2(wi,0, 

which turns into ([t]) in the hmit t 10. See Section [sTT] 

• The scheme is appHcable to arbitrary initial data u'^' G L^{I) with finite energy. In fact, 
weak convergence of to vP suffices to conclude strong convergence of the discrete 
solution MA in i^([0,T] x /) to the correct weak solution u of ([T]). See Theorem [l] below. 

• Discrete solutions obey a minimum/maximum principle. See Section [5.2[ 

The choice of D^^(/) originates from an alternative formulation of equation ([T]). Namely, u is a 
(positive and classical) solution to ([T]) iff its inverse distribution function X — see Section [2] for 
its definition — satisfies the initial-boundary value problem 

atX = V'(Xe)«-y.oX, X(i;0) = 0, X(i;M) = l, X(0; = X"(C), (10) 

where : M+ ^ M is defined by 

V;(s) = s0(s-i) for all s > 0. (11) 



The variational structure of (10) is quite apparent: solutions X to (10) are gradient flows of the 
functional 



E(X)= V(X5(0)d$+ / VoXiOd^ (12) 
Jo Jo 

with respect to the usual scalar product on L^([0, M]). In effect, we discretize the L^-gradient flow 
(10) rather than the W2-gradient flow ([T]), representing X as a linear combination of piecewise 
linear ansatz functions with respect to the (time-independent) mesh ^. 

1.3. Convergence result. Our main result is the following. 

Theorem 1. Let a non-negative initial condition u'^ G L^{I) of mass M with E(u*') < oo be given, 
and fix a time horizont T > 0. 

Consider a sequence o J discretizations A'^-'' = {t^''^;$}-^^), consisting oj time steps r^^' \. and 
spatial meshes ^^■'^ with maxfc((5^"''') \. 0, and an associated sequence of initial conditions u^^fj^ G 
D^^)(J). Assume that u^^^j J — >■ m° weakly in (I) , i/iat E(-u^(j) ) < E, that maxk 5[^^ / mini S'f^ < 
a, and that the following inverse CFL condition holds: 

(max^^'V < eV-'Y '^"''T V^'^ (13) 

with tjj defined in (11 ), and with A > from ([6|. 

The scheme ([9| produces a sequence of discrete solutions u^u) . Denote by u/^u) '■ [0, oo) — ^ 
Dt^(/) the respective interpolants that are piecewise constant in time, see (62). Then u^u) con- 

Sj 1 1 

verges strongly in Li([0,r] x /) to the unique weak solution u of 



A comment is due on condition (13). Since ip'^s) — for s — > oo, this condition implies that 
the non-negative initial datum needs to be approximated by strictly positive data , and 
the smaller one whishes to choose the minimal value of , the finer one needs to make the grid 
Condition (13) thus quantifies the intuitive requirement that not only the mesh of S^^s in 



[0, A/], but also the induced mesh of Xfc's in / should become arbitrarily fine in the limit, uniformly 
for all times t G [0, T]. Consequently, our scheme does not allow to track propagating fronts — 
like spreading Barenblatt profiles — directly on the discrete level as in O [T^ , but it is able to 
approximate these fronts arbitrarily well with strictly positive solutions if the mesh is sufficiently 
fine. 
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1.4. Related results from the literature. Studies on Lagrangian schemes for ([ij are widely 
scattered in the hterature. Aheady MacCamy and Sokolovsky [I6j present a discretization that 
is almost identical to ours, for (IT]) with P(u) = and V = 0. Another pioneering work in this 
direction is the paper by Russo [21], who compares several (semi-)Lagrangian discretizations in the 
linear case P(u) = u; extensions to two spatial dimensions are also discussed. Later, Budd et al [5] 
used a moving mesh to capture self-similar solutions of the porous medium equation on the whole 
line. The general theme was picked up recently by Carrillo and Moll [7], who define a Lagrangian 
discretization of aggregation equations in two space dimensions, based on the reformulation in 
terms of evolving diffeomorphisms [lOj . 

The connection between Lagrangian schemes and the gradient flow structure of equation (II]) 
was investigated by Kinderlehrer and Walkington TJ] and in a series of unpublished theses [111 Ho] . 
In a recent paper by Westdickenberg and Wilkening [21| , a similar scheme for (|T]) is obtained as a 
by-product in the process of designing a structure preserving discretization for the Euler equations. 
Burger et al [H] devise a numerical scheme for (JlJ in dimension two on basis of the gradient flow 
structure, using the hydrodynamical formulation of the Wasserstein distance [3] instead of the 
Lagrangian approach. The Lagrangian approach was adapted to fourth order equations, namely 
by Cavalli and Naldi [5] for the Hele-Shaw flow, and by Diiring et al [5] for the DLSS equation. 

In the aforementioned works, numerical schemes are deflned and used in experiments; qualitative 
properties and convergence are not studied analytically. Some analytical investigations have been 
carried out by Gosse and Toscani [12] : for a Lagrangian scheme with explicit time discretization, 
they prove comparison principles, and they rigorously discuss stability and consistency. Also, a 
full discretization of the Keller-Segel model has been analyzed by Blanchet et al [3] in view of 
convergence to equilibrium. However, to the best of our knowledge, a proof for convergence of 
discrete to continuous (weak) solutions is not available in the literature. 

Finally, a remark is due on an alternative way of proving convergence of the scheme. By use 
of stability results for gradient flows [TJ [2] and the machinery of F-convergence, it seems likely 
that Theorem [T] can be obtained by exploiting the variational structure more deeply than we do 
here. In particular, the theory on perturbed A-contractive gradient flows developed by Serfaty 
[22] indicates an alternative route towards the same goal. We followed the elementary approach 
based on a priori estimates here, partly in order to avoid heavy machinery, but mainly with the 
aim to develop a "stable" concept of proof that generalizes more directly to gradient flows without 
convexity properties (like fourth order equations). 

1.5. Outline of the paper. Section[2]below summarizes some basic results on inverse distribution 
functions and convexity in the Wasserstein metric. In Section [3] we describe in detail the spatial 
discretization and study the restrictions of the Wasserstein metric and energy to the D^^ (/) . The 
discrete scheme (|9| is studied in Section |4] and we derive the Euler-Lagrange equations. Section 
[5] provides a summary of some qualitative properties of the discretization, like metric contraction 
and the minimum/maximum principle. The proof of Theorem Jl] is given in Section |6] The paper 
concludes with the results of various numerical experiments in Section [7] and with a calculation 
of the consistency order. 

2. Preliminaries and notations 

For an introduction to the theory of optimal transportation, we refer to f23J . A comprehensive 
theory of gradient flows in the Wasserstein metric can be found in [Ij. 

2.1. Inverse distribution functions. Throughout the paper, we shall denote by 
D*^(/) := |ii e L\l) n L°"{I) ess^nf w(a;) > 0, ju{x)Ax = A/| 

the space of positive density functions of total mass M. For u E D^^{I), define its distribution 
function U : I ->-[0,M] by 

U{t;x) = / u{t;y)dy, 
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and introduce its inverse function X ~ U ^ : [0,M] /. By our choice of D*^(/), the latter is 
weU-defined and belongs to 

X {X e C°^^([0, M];I) \ X(0) = a, X(M) = b, X strictly increasing}. 

Thanks to the Lipschitz continuity of U and X, we can differentiate the identity U o X(^) = ^ at 
almost every ^ e [0, M] and obtain the relation 

u(X(e))X^(0 = 1 for a.e. ^ e [0, M]. (14) 

The inverse distribution function allows for an explicit representation of the Wasserstein distance 
in one spatial dimension. 

Lemma 2 (see e.g. [53]). Let uq,ui E D*^(/) have inverse distribution functions Xo,Xi e X. 
Then their Wasserstein distance amounts to 

f rM \ 1/2 

w2(uo,ui)- / [xm-xo{o?da , (15) 



and a minimal geodesic (us)o<s<i connecting uq to ui in W2 is given by Xs = sXi + (1 — s)Xq. 

2.2. Properties of the energy. Given P : IR>o — ?> M, let : M>o -> M be an arbitrary second 
anti-derivative of r P'{r)/r, and define ■0 ■ ^+ — >■ R by ip{s) = s0(l/s). Introduce the 
functionals E on D^^(/) and E on X, respectively, by ([s]) and ( [T2| . 

Lemma 3. For every u G D*^(/) with inverse distribution function X € X, one has 

E(it)=E(X). (16) 

Further, (j) is strictly convex and satisfies 

P(r) = r(j)'{r) + 0(0) - 0(r) for all r > 0. (17) 
Finally, ■(/'(•s) — >■ 00 for s ], 0, and ip" is a positive non-increasing function. 

Proof. We perform the change of variables x — X(^) under the integrals in the definition (|5| and 
use ( [14| : 

E(u)=/ 0(w(x(O))X4(Ode+ / i/(x(OMx(e))x^(e)d^ 

Jo JQ 

(•M , T , rM 



I Kx^y^^^^'^^^ I nX(0)de = ]E(X). 



The claims about (j) and V' are direct consequences of the hypotheses in ([2]). By definition, (^"{r) — 
P'{r)/r > for all r > 0, so cj) is strictly convex. ( 17) follows by differentiation of both sides w.r.t. 
r > 0. It follows further that 

^'(,s) = 0(s-i) - s- V'(s-') = <^(0) - P{s-'), 
^"{s) = s-^P\s-^) > 0, 
V;"'(s) = -d2p(s-l)/ds2 >0, 
SO ip" is indeed positive and non-increasing. Finally, 

limV'(s) = V(l) + lim / P(cr"^)dcr = -0(1) + lim / dp = +00 

slQ slO r^cc 

since P'(p) — > 00 for p — !• 00, and hence also P{p)/p — > 00. □ 

The convexity of the functional E with respect to the Wasserstein metric is most conveniently 
studied when the latter is considered as a functional of X instead of u. Indeed, by Lemma [2] 
above, geodesic interpolation between uo,ui € D*^(/) corresponds to linear linterpolation between 
Xg, Xi G X. 
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Xi X2 X3 




^1 ^2^3 



Figure 1. A typical density function u G D^{I) (left) and inverse distribution 
function X G X. 



Lemma 4. The functional E is bounded from below, 

M 



E(X) >E:= ib-a)4)(—-) +MinmV(x), 
and it is {—A) -convex on X with the A given in i.e., 



for all X'^jX"'^ e X, and every s G [0, 1]. 

Proof. Since is convex, the lower bound follows by Jensen's inequality: 



(18) 



E((l-s)X" + sXi) < (l-s)E(XO) + sE(Xi) + — ^ ^/ [X°(0 - Xi(e)]2 df (19) 



By definition of V', this yields ( 18 ). Next, let X°, X G X and s G [0, 1] be given. Since tp : 



is convex by hypothesis, it follows in particular that 

^{1 - s)x^(c) + sx\[o) < (1 - / ^{^{0) de + W V'(x|(e)) de 



Further, a Taylor expansion yields 



A 



V{{1 - s)y + sz) < (1 - s)V{y) + sV{z) + -s(l - s){y - zf 



for arbitrary y,z I. In combination, this implies inequality (19 1 



□ 



3. Spatial discretization 

Inside the space X of inverse distribution functions, we define the finite-dimensional subspace 
X| of those functions, which are piecewise affine with respect to a given partition ^ of [0, M] 
into sub-intervals. Correspondingly, there is a finite-dimensional submanifold D^^(/) of D*^(/) 
consisting of those densities, whose inverse distribution functions belong to X^. Densities in D|^(/) 
are piecewise constant. Since we shall work simultaneously in the spaces D|^(/) and X|, we need 
to introduce various notations. 
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3.1. Ansatz spaces. A vector ^ — {£,0,^1, ■ ■ ■ ,£k) with entries £j such that 

= $0 < a < • ■ • < = M 
defines a partition of [0, M] into K sub-intervals. We denote the lengths of the intervals by 

4 = £k - S.k-1 for aU fc = 1, . . . , X, 

and introduce further 

^(0=min4, J(0=niax4, «(0 = ^fv- (20) 

k k d[^) 

The associated {K — l)-dimensional reduction of the space X is given by 

:— {X e X I X piecewise affine on each [Cfc-i, for A: = 1, ... , X}. 

Functions in are conveniently represented as linear combinations of the K +1 hat functions Oq 
to 9k defined by 

{{S. - S.,n-l)/S,n for < C < Cm (if TO > 1), 

(6«+i - O/^ni+i for < C < Cm+1 (if m < if - 1), 
otherwise. 

More precisely, there is a one-to-one correspondence between functions X G and vectors in 

y := {x = (xi, . . . , xk-i) I a < xi < a;2 < • • • < xk-i < b} C I^^^; 

this correspondence is established by means of : y — >■ X| , with 

K 

X = X^[x] = J2xk0k- (21) 

k=0 



Remark 5. By definition, x e y has components xi to xk-i- In (21 1 above and for the rest of 
the paper, we shall always use the convention that 

Xq — a and xk = b. (22) 

Also, we introduce in analogy to (5(^) the mesh width 

5(x) max(a::fe - Xk^i). (23) 

k 

Occasionally, it will be more convenient to work with vectors z G M.^ of difference quotients: 
define : y by 

z = (zi,...,z;f) = zg[x] gM^, with zfc^^^i— (24) 



using again our convention (22 1. 

Finally, we introduce the associated (A' — l)-dimensional submanifold D|^(/) := u^[y] C D*^(/) 
as the image of the injective map : y — ^ (/) with 

= where = (25) 

k=i ^'^-i 

3.2. Representation of the Wasserstein distance. The Wasserstein distance between any 
two elements of D^(J) is easy to compute using (151. 



Lemma 6. Fix a discretization ^, and let u ,u G D| (/) have representations u = Ug[x^], 
= u^[x-'^], with ^,x^ G y. Then 

W2(m",u1)2 = (xO-xi)^W(x°-xi) (26) 
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with the symmetric tridiagonal matrix W = {y^m,k)m k=i ^ i) given by 



62 2{52+53) 





V 

Moreover, for every v £ M.^^^ 

K-l 



Sk-1 2(6k-i+6k)J 



(27) 



A'-l 



^ 51 ('^'^ + ^k+i)vl < v'^Wv < ^ ^ (4 + 5k+i)v 



(28) 



k=l 



fc=l 



Proof. Plugging the representation (21) into the definition (151 of the Wasserstein distance yields 



Since Xq = Xq — a and x'jf = x]^ — b, the last sum actually runs only over indices from one to 
K — 1. Therefore, in order to prove (26), it suffices to show that 



M 



for all m, fc — 1,...,K — 1. Since dm has support [(,m-i, £.m+i] by definition, it follows that 
Wm fc = if Ito — /cl > 2. Moreover, for m = fc we have 



M 



and for fc — m + 1 , 

rM 







1]'^ dr] + Srn + l C^dC = -{Sm+Srn+l), 



Finally, let v E M^^^^ be given and observe that 



(fm+1 - 0(C - 6n) dC = Sm / (1 - V)'n drj = -Sr, 



K-l 



K 



m=2 



ni—1 



From here, (28 1 is immediately deduced using binomial formulas. 



□ 



3.3. Representation of the energy. The restriction of the energy E from (16) to the subspace 
X| is naturally associated to the functional : y ^ M with 

Eg(x) :=E(X^[x])-E(u^[x]). 

By straight-forward calculations, one obtains the following more explicit representation. 

Lemma 7. For every x e y, we have 



K 



E^(x)-^4V^ 



fc=i 



Xk - Xk-l 



F(x^[x](0)de. 



(29) 
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Moreover, the (Euclidean) gradient vector dsiE^{x) = (^dx^^^^{x)) 



[9xE^(^)], 



Om+1 



K-1 
m—1 
M 



pK-1 



is given by 



id the Hessian matrix 9|E|(x) — (^d^^xk^^i^)) 



K-l 
m,k—l 



'(K~i)x{K-i) symmetric with 



(30) 



-r 



rn.k 



t5m + l 



M 



(31) 



M 



Vxx{'X-iM)0rn0m-l if k ^ m - 1, 

if 1 < k < m ~ 1. 



Further, the functional inherits boundedness and convexity from E. 



Lemma 8. E^ is bounded from below by'E defined in (18 1. Further, it is {—A) -convex with respect 
to the quadratic structure induced by W , i.e., V'^E^(x) + AW is positive semi-definite for arbitrary 
x° e y. Consequently, 

(xi - xO)^(a5cE|(xi) - 9jE^(xO)) > -A(xi - x°)^W(xi - x°) (32) 
holds for every x°,x^ € p. 



Proof. Boundedness from below is a trivial consequence of ( 18 1 and the definition of E^ by re- 
striction of E. Convexity is a direct consequence of the convexity ( 19 ) of E, taking into account 
(26), and that is an affine map. The estimate (32) is obtained by Taylor expansion. □ 

4. Time-discrete evolution 

Throughout this section, we fix a pair A — (r, ^) of a time step with r > and a spatial 
discretization ^ = (^g, . . . , S,k)- 

4.1. Minimizing movements. With the finite-dimensional manifold D|^(/) given at the end of 
Section 3.1 above, the procedure ^ can now be used to define inductively — starting from a 
prescribed initial datum e D^^(/) — a discrete solution ua ■= ('UA)n^o- 

Proposition 9. Assume that tA < 1, with A > defined in ([6|. Recall the definition ofEr from 
([S]). Then, for every G D|^(/), there is a sequence {u^)^^q, such that u\ € D|^(/) is the 
unique minimizer of 'Eir{' ,u^^) on the restricted set D|'^(/), for every n G N. 
Moreover, define the associated sequence (x^)^o ''/^ ^ f by 

u^[x^]. (33) 



-'A 



Then each x^ is the unique solution x G y the system Euler- Lagrange equations 

1 



:w(x-x'A-^) 



(34) 



with c)5jE^(x) explicitly given in ( |30[ ). 

Proof. By definition of D|^(/) as the image of p under u^, it suffices to prove unique solvability 
of the minimization problems 

Ea(x,x^-1) := ;i(x-x^-i)^W(x-x'A-i)+EA(5) ^ min 
It 

for X G p. To this end, observe that 

EA(f,2'A"^) =E€(x) + |(f-^-^)^W(x-xA-i) + i(r-i - A)(x-xl-i)^W(x-xA-i) 

for every x G p. From Lemma [4j we know that the sum of the first two terms on the right-hand 
side constitutes a convex function in x G p. Since rA < 1, and since W is positive definite by 
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LemmajBj the last term is strictly convex. Thus, EaI'jX^^ ) possesses at most one critical point 
in p. 

To show the existence of a minimizer, let (x'-'')jgN be a minimizing sequence for Ea(-,x"^^) in 
y^. Since each of the K — I components a: belongs to the compact interval /, we may assume 
without loss of generality that x*^^^ converges to some x* G I^^^. It remains to be proven that 
X* € pg. Since (x'^''))jgN is a minimizing sequence, Ea(x*^-'-',x"^^^) is bounded, and so, for every 
me {!,..., if}: 



C>l(xO) 



)^W(x 



U) 



-<n-l 

X A 



K 
k=l 



M - S„ 



4 



M min^fa:;) 

xei 



M 



y(x«[xW](0)dC 



Since ^pis) — > oo for s 4- 0; this implies that — x^_i > eSm > with some e > for all n G N, 
and thus also x^^ — x*^_i > eSm > 0, implying x* G y. By continuity of Ea(-, x27^) on y, it follows 
that X* is a minimizer. 

The argument shows that Ea(-,Xa^^) possesses a unique critical point in y, thus the corre- 



sponding Euler-Lagrange equations (|34|) are uniquely solvable. □ 

4.2. Euler-Lagrange equations for the difference quotients. The analysis that follows will 
be based primarily on another representation of the system (34 1 of Euler-Lagrange equations, 
which is formulated in terms of the difference quotients = (z", . . . , z^) introduced in (24 1: 



(35) 



To begin with, we introduce quadratic analogues 8i, . . . , Qk ■ [0, M] — > M of the piecewise linear 
hat functions . . . , 9x as follows. Let the numbers 71, ... , 7^^ G ( — 1,1) be defined by 71 = — 
0, and 



7fc 



5k+i — 5k-i 
4+1 + 24 + 4- 



for k 



,K-l. 



(36) 



Then the Qk are given by 



efc(0 



for fc = 2, 



= < 



^(4+1 + 4 + 4-1) - 4^(2^ - (^fc + ^^^-i) - 7fc4)^ 



) 

if — 1, and by 



if a-2 <^ <a-i, 

if^fc-i <e<a, 

if C/c < C < &+1, 
otherwise 



ei(e) 



l_f2 



25i 







2(5 



ifO<e<Ci, 

if6<e<6, 

otherwise, 







otherwise. 

Lemma 10. For each k — 2, . . . , K — 1, the function is supported on [S,k-2y£,k+i] o-nd satisfies 

-{Qk)i = {i-ik)ek~{i + ik)ek-u (37) 

and we have (61)5 = ~0i and (6/^)5 ~ &k-i- 



Proof. This follows directly from the definition. 



□ 
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Next, we define the matrix W = O^m.m') 



K 



t,KxK 



by 



(38) 



The matrix W essentially plays the same role for the as W for the x^. Its entries are more 
complicated, but still can be calculated explicitly. 

Lemma 11. The matrix W is tri-diagonal and has entries 

^ ¥1, + ^5^5^+1 + i^<5™<5™_i if2<k^m<K~l, 
1a2_l1a.a. ifk = m = l, 

if k ^ K, 
if k = m — 1, 



< 



^Sf + ^6261 

6 "m+1 





if k = m + 1, 
otherwise. 



(39) 



Proof. The explicit representation of W is obtained by a tedious, but straight-forward computation 
that requires nothing but integration of quadratic polynomials and the use of ( 36 1 . □ 

Lemma 12. For the solution (x^)5^o obtained in Proposition^ let (z^)5^o associated 
sequence from ( |35[ ) . Then each satisfies the following system of Euler-Lagrange equations: 

1 

T 



[W(Z-A - ^r')]ra = (1 - 7™)^'(C+l) - 2V'(C) + (1 + 7™)V''(C-l) 

m+1 ^Cfc (40) 



^ / K4xg[x"])e„,de 



for every m = 2, 



, -fC — 1, and 



-[w(zA - zr')]i = ^'(4^) - V''(^r) - E r v,,{x^[^])e, d^, 

''" fc=i •'ik-i 

J[w(zA - zX')]a' = i^'izK-i) - ^'(zk) - E K.(X«[x"])eKdC. 



(41) 



(42) 



Proof. Fix an index m S {2,...,_fC — 1}. With 7^ given by (36), multiply the mth and the 
(ra — l)th component of the Euler-Lagrange system (34) by 1 — 7™ and 1 -|- 7„i, respectively, and 
substract the latter from the first. This yields 



1 /l-7„ 



1 - 7™ c 



(Om + Om+i) 0„ 



3 

1+7. 



m— 1 
1 1 ^ri 



= (1 - 7™)^'(4\+i) - 2V/«J + (1 + 7™)V''«^.i) 

pM 

- / K(X^[X"]) [(1 - 7„0e™ - (1 + 7m)^m-l 

Jo 

The expression on the left-hand side can be rewritten as 
1-7 



l+7»> 



de 



3m— 12;,„_2 



6 



1 - 7m . 1 - 37„ 
Om+1 H 7. On 



1 + 7m . 1 + 37m 
?i Om-l H On 

2 6 



6 

l+7r- 



''m-l 



6 



"<5m-i(2;„_i 2;„j_2) 
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where we have used the relation 



1 - 7m . 1 - ilm f. _ 1 + 7™ . 1 + 37,„ . _ 1 ^r. 

n "m+1 i Om — n Om-1 i — -^^ — y\m,m, 

2 6 2 5 dm 



which is a consequence of our definition of 7^ in (36 1. Thus, we obtain 
1, 



:[W(zl-zl 



= (1 - 7m)V''(C+l) - ^i^'i^ln) + (1 + 7m)V''(2™-l) 
/ 

- / 14(X5[xl])[(l-7„)0„-(l + 7™)0™_i]de. 
Jo 

Using property (371 of the function 0„i, and integrating by parts, we arrive at 

- / K(X^[^A]) [(1 - 7™)^m - (1 + 7m)^m-l] = / V,{y.^[l\]){Q,^)^ d^ 

Jo Jo 



(43) 



M 



(44) 



(XcK]),^..(x^[xi])e„,de 



The boundary terms in the second hue vanish, since 9^(0) — 0m(M) — 0. FinaUy, observe that 
(X^[x'^])^(^) = Zk for aU < ^ < ^fc. Insert the result into ^ to arrive at 



Equations (|4l|) and (|42|) are directly obtained from (|34|) for m ~ 1 and for m = K — 1, 

of 

□ 



respectively, after an integration by parts like in (44 1. The boundary term vanishes because of 
hypothesis (|3]). 



4.3. Energy dissipation. The estimates derived below are at the core of our convergence proof 
in Section[6j We start with two energy-type estimates, which are classical in the theory of gradient 
flows. We recall that u\ ~ u^[x^], and that the gradient 9xEg is explicitly given in (30). 



Lemma 13. For every G N, 



1 ^ 



n=l 



N 



(45) 
(46) 



Proof. Since minimizes Ea(-,uA ), we have in particular Ea(maj'"a ) — Ea(uA ) 
which implies that 



2t 



W2«, <-'r < E«-0 - E«) 



Evaluation of the telescopic sum yields ( 45 ) . To obtain ( |46| , multiply the system ( 34 ) of Euler- 
Lagrange equations by (t/2)^/^W~^/^, take the Euclidean norm on both sides, and sum over n = 1 
to n = TV: 

T ^ ^ / 1 



Insert this in (45) to obtain (46). 



□ 



Proposition 14. For every N £N, 

Sk + 4+1 



E E ^"^^^^L^^^^ ^ "(^) [E«)-E«) + MTsup(K(x)^)], (47) 



n=l fc=l 



xei 



with T — Nt and the ratio a(^) being defined in (20). 
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Proof. From the upper estimate on W in ( 28 1 , it follows that 



> m-' E 



fe=i 



M 



1 2 



1 



fe=l 

where we used that {x — j/)^ > a;^/2 — for arbitrary a:, y g M. Now, on one hand, 
E / K.(X^r])0fedd <sup(\4(a;)') E / ^'^'i^ 



M \ 2- 



k=l 



And on the other hand 

K-l 



k=l 
K-l 



(4 + 4+i)' 



k=l 



<M<5(0sup(K(x)^ 



Combining these estimates, we obtain 

E - f ^[5x^%(^1)]"W-M9.E,(x1)] + Msup (VAX) 



KG/ 



Multiply by r and sum over n = 1 to n = A^. An application of (46) yields (47). 

5. Qualitative properties of the discretization 



□ 



Throughout this section, we fix a space-time discretization A = (r, ^) and consider a given 
discrete solution ua = {ua)^=o- 

5.1. Metric contraction. One of the fundamental properties of our solution scheme is the preser- 
vation of the contraction property ([7| . 

Proposition 15. If va = {va)'^=o '^"■2/ other discrete solution, then 

W2{ul,vir < (l-2Ar)-«W2«,0' (48) 

for all n eN. 



Remark 16. Since (1 — 2At)" < exp(— 2Anr) for every n G N, estimate (48 1 is slightly worse 
for every t > than the limiting estimate ([T]). 

Proof. For xa, Ya such that Ua = "^[x^], v^ = u^Iya] we know by Proposition |9] that 

W(xl - x'^-i) - -ra^E^l^l) and W(yl - y^') = -rds^E^irA) ■ 
Substracting these equations, we obtain 

Wi/2(x'A - yl) + rW-i/2(a.E^(x^) _ d^E^{rA)) = W^/^^^-i _ ^^-i)^ 

where W^/^ and W~^/^ are the (symmetric and positive definite) square roots of W and its inverse, 
respectively; see Lemma [6j Taking the norm on both sides yields 

(x'A - y2)W(x'A - y^) + 2r(xA - yA)^(9,jE^(xl) - ds^E^if^)) 



A 



Combining this with the convexity property ( 32 1 , we arrive at the recursive relation 



(1 - 2At)(xi - yi)^w(x'A - yi) < (^r' " frv^i^r' - rr')- 



Iteration of this estimate and application of ( 26 ) yields ( 48 ) . 



□ 



14 



DANIEL MATTHES AND HORST OSBERGER 



5.2. The maximum and minimum principles. Recall that ■0" is a positive and non-increasing 
function by Lemma [sj and recall the definition of a(^) from (20 1. 

Proposition 17 (Minimum principle). Assume that At < 1/2, and assume further that ^ and t 
are related by the inverse CFL condition 



d{S,)^ < 6ip" (Z^)t, where Z^, := 



6a(^)e 



miUj. u 



2TA 



Then, for every n with nr < T , 



minu^ > e 



-2nTA , 



(49) 



(50) 



Proof. The minimum principle ( 50 1 for ua is equivalent to a maximum principle for za . By 
induction on n, we prove 



maxz^ < (1 - Ar)"" maxz^, 



(51) 



which is a slightly sharper estimate than (50) since (1 — rA) > e^^^^ under our assumption 
tA < 1/2. So fix n with ut < T and assume that < Z^""^) for all k. Suppose that 

Z'"^ = for some m € {2, . . . , if — 1}. We estimate the integral term in the Euler-Lagrange 
equation (|40]), using the positivity of the and that 9m > 0: 



[Wz"]™ < [Wz'^-'U + r((l - 7™)V''(C+i) - 20' (O + (1 + 7m)^'(C_i)) 

m+l 



At ^ z^ / e„,de 



Recalling (38), we conclude further that 
(1-At)[Wz"],„ < [Wz"-i]„, 

+ (1 - 7m)r(0'(z;;+i) - i^'izl^)) + (1 + 7™)r(V'«n-i) " V''(0)- 
Since Z^") = > for all fc, and since -0" > is non-increasing, it follows that 

^'(c+i) - ^''(o < -0^"(z("))(z:j, - c+i), 

V^'(C-i) - i^'i^l) < -0"(z("))(C - C-i)- 

In particular, these terms are non-positive, and so (52) implies 

(1 - At)W,„,™Z(") < (T,„Z("-i), with a,„ = W,„,„,_i + W„,,„ -f W 



(52) 



(53) 



m,m+l ■ 



From the explicit form of W in (39), we obtain am < 3a(^)Wm,m, and by means of the induction 
hypotheses, we conclude the rough bound 



< 6a(|)Z"-i < 6a(|)(l - At)'^"-^^ maxz° < Z, 

k 



T- 



(54) 



We return to (52), insert (53), use that ip" is non-increasing, and find after some manipulations: 

(1 - AT)a™Z(") < amZ("-i) + (1 - 7™)(^^e+i - ^^"(^J)) " C-i) 

1 - At 



+ (1 + 7 m) 



6 



The inverse CFL condition \4Q\ implies non-positivity of the last two terms, so (|54| refines to 



^(") < (1 _ At)-1z("-i) < (1 - At)-" max 2? 



If the maximum is attained at one of the boundary points, m — 1 oi m = K, then a similar 
calculation can be carried out using ( 41 ) or ( 42 ) , respectively, instead of ( 40 ) . □ 

Similarly to the minimum princple, one obtains the following maximum principle. Notice that 
the inverse CFL condition is the same as before, i.e., it involves min^. and not max^. w^. 
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Proposition 18 (Maximum principle). Let all the hypotheses of Proposition 17 hold, and assume 
— in addition to (49) — that 

(1 + At)^(^)2 < 6tV'"(Z;), where A:=maxK^>0. 



Then, for every n with nr < T, 



maxw^ < e^"^'^maxu^. 



We omit the proof, which is very similar to the one given above. 
5.3. Regularity. Introduce the function $ : M+ — > M by 

^r)^ f ^^"{p)dp. (55) 



$ is strictly increasing, with 

$'(r) = v^0"(r) = r^^^^ij"{r-^) (56) 



by definition of V' in (11). Further, $ grows asymptotically faster than P, 

hm -y- = +(X), (57) 
which is easily verified by I'Hospitals rule: 

lim -VV = 2 lim ^ , y = 2 lim —-d- = 4 lim (Vf$'(r)) = 4 lim P'(r) = +oo, 

r—^oo P(7') r— J-oo P'(r) r— >oo WT' r— fcxD r— fcxD 



due to ([2]). In particular, $^ is superlinear at infinity. 
Proposition 19. For every N ^N, 



N 

r5]Var^$(^")2) <c(E«(^o^),a(0,A^T), (58) 

n=l 

where 

/ M s 2 

C(£:,a,r) = T$( ) +6(b^ a)an+a)\E + MTsup{VJxf)]. (59) 



Estimate ( 58 ) mimicks the standard energy dissipation estimate for (jlj , which is 
-^E(7.) = i.[0'(7.) + V]ldx> ^u)l dx - C. 
Roughly speaking, the integral of the spatial derivative is replaced by the total variation, 

A'-l 

Var^ {^{u\f) = l*K+i)' - 

k=l 

which is about the maximal regularity that can be associated to a piecewise constant function. 



Proof. Recall that — 1/z^ = 6k /{x^ — x'^_^). On one hand, using (56), 



2 



($«+i)-ci>«))^= (^1^^ a>'(r)dr^ 

= (^^"'^'r-5/2^"(r-i)dr^ =(^^'^"' ^s^" {3) d^ 

< max(z,", 4Vi) (£' r{s) ds^ = max(z,", Zfe"+i) (^'(z,") - ^'(4Vi))' 



On the other hand, recalling the definition of a(^) in ( |20[ ), 

(4 + 4+i)(4 + 4+1) < (1 + «(0)(44 + ^Vi-^fc+i)- 
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We combine these estimates and sum over k — 1,...,K— Ito obtain 



Since xJJ^i ~ = 2^'(5fc + z^'_|_i'^fc+i, it follows further that for arbitrary fc, ^ e {1,2,..., K} with 
fc < 

£-1 

|$«) - $«)| < ^ - 

m— 

< (2(6 -a)5")'^'. 

Since u has average value M/{h — a), there exists an index m* with Um* < M/{b — a), which 
implies the absolute bound 

max <^{uk) < <i>K„0 + - <f{ur^-^)\ < <^(j^) + (2(6 - a)S"Y^\ 
i<k<K \b — a/ 

So, finally, 

K-l K-l 

Var^ {<P{uir) = ^ |<i>«+i)2 - $(0^1 < 2^max^<i>«) ^ - $(«^)| 

fe=i ^ ^ fe=i 



< 2 



(^) + (2(6-a)5")^/^](2(6-a)5")^ 



M \2 

<*(l ) +6(6-a)5. 

^ — a/ 



To obtain (58), sum with respect to n = 1, . . . , and use the energy estimate (47 1. □ 

6. Convergence 

In this section, we prove Theorem [ij Let a time horizont T > and an initial condition u'^ e 
L^{I) with E(w°) < cx) be given. We consider a family Aj = [Tj^^j) of time-space discretizations, 
with j e N. Accordingly, we denote by Kj the number of nodes of and Nj is the smallest 
integer with tjNj > T. 

Throughout this section, we assume all the hypotheses of Theorem [TJ 

• Tj 4, and i as j — 7> oo; 

• initial conditions u^. e D|^(/) are given for each j, such that vP^, converges to vP weakly 

inii(/). 

• uniformly in j G N, 

— - / fi?>p2Ar \ 

a(l,)<a<oo, E«^.)<E<oo, <5(^^f < 6^" ^ r,. (61) 

\ mm^; J 



Denote by (wAj)n^o the corresponding discrete solutions obtained as in Proposition |9j and intro- 
duce the time-interpolated functions : [0,T] — D^(/) by 

"A, {t; x) = u^^. (a;) for all te{{n- 1)tj, nr^] n [0, T]. (62) 
The following preliminary result plays an important role in the convergence proof. 



Lemma 20. With the maximal mesh width (5(x) defined in (23), we have 



max (5(x^ . ) — > as j —> oo. (63) 
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Proof. This is a consequence of the minimum principle and hypothesis (61 1. It follows from the 
assumption limr^oP'(^) < oo in ([2]) that 



o2AT 



mm^ u\ . 



for some appropriate constant C. Recalling that (xk — Xk-i)ui; = 6k, Proposition 17 thus implies 
— uniformly in j and n < Nj — that 



J \ TTHiTi I \ 1-mn i/V / 7/;'' ( ^ 1 

^ V min^ u*! . / 



using hypothesis (61 ). Since Tj ,|, as j — >■ cxd, this proves the claim. □ 

6.1. Compactness in W2. The following weak convergence result is a well-known consequence 
of the energy estimate (45 1 in combination with the Arzela-Ascoli theorem. 

Proposition 21. More precisely, every subsequence of {uAj)jeK contains a sub-subsequence that 
converges uniformly w.r.t. t G [0,T] in W2 to a limit curve £ C"'^/^([0, T]; W2). 

A proof can be obtained by application of ,1,, Proposition 3.3.1]. 

6.2. Compactness in L^. The following compactness property on the ua^ is at the basis for our 
convergence proof. 

Proposition 22. Every subsequence o/(uAj )jeN contains a sub-subsequence such that the respec- 
tive UAj converge to some u*, and the P(uAj) converge to P(u*), both strongly in L^([0,T] x /). 

The proof of this proposition is an application of the Aubin-Lions compactness principle. Specif- 
ically, we use: 

Theorem 23. [Adapted from Theorem 2 in [50]/ Assume that: 

(1) There is a normal coercive integrand ^ : L^{I) — [0, 00], i.e., 5" is measurable, lower 
semi- continuous and has compact sublevels in L^{I), for which the following is true: 

sup / 5(uA, (t))dt < 00. (64) 

JGN Jo 

(2) The UAj are integral equicontinuous with respect to W2, 

limsup / W2{uA,{t + h),UA^{t))dt^Q. (65) 

hit) j gpj Jq 

Then the sequence {uAj)j&i "i-s relatively compact in L^{[0,T] x /). 

In order to define ^, we first recall that the total variation of a function / G L^{I) is given by 



{/}tv '^'^P I / ■ 



f{x)Lp'{x)dx 



^GCi(/),sup|^(x)| <1 



xei 



It is easily checked that for piecewise constant densities u e (/), and with $ from (55), we 
have 



Ki-l 



{$(u)2}tv = Var^ ($(w)2) = |$(«,+i)2 - <^{ukf\. (66) 



k=l 



Now let 5^ : L^/) -> M U {+00} be given by 

f{$(u)2}Tv ifwGD*^, 



-00 otherwise; 



where D*^(/) denotes the closure of D^'^(/) in L^{I), which consists of all non-negative L^-functions 
with integral equal to M. 
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Lemma 24. The functional J? defined above is lower semi- continuous and has relatively compact 
sublevels. 

Proof of Lemma^ Let A, := ^'-^((-oo; c]) C L^{I) be a sublevel of S^. By HH Theorem 1.19], 
the set Be := | u G Ac} is relatively compact in L^{I); here we use that our domain / is an 

interval, so that {$(u)^}tv < c and u(a;) da; = M induce a uniform bound on the BV-norm of 

Thus, if {ui) is a sequence in converging to uq in L^{I), then also converges to 

$(wo)^ in L^{I). By lower semi-continuity of the total variation {-Itv [HI Theorem 1.9], the lower 
semi-continuity of ^ follows. 

To conclude compactness of Ac, it suffices to prove that the mapping u is L^{I)- 

continuously invertible. For that, let a sequence (/^)fgN in be given, which converges to 
some /o in L^{T). Since the map r i— ^ $(r)^ is strictly increasing, positive, and continuous with 
superlinear growth (57), it possesses a strictly increasing, positive and continuous inverse with 



sublinear growth. Hence, there are a uniquely determined sequence of functions ui € Ac such that 
^{uiY = fi for all £ e N, and a unique uq G Ac with <i>(uo)^ — fo- We wish to show that ug 
converges to Uq in L^{I). By standard arguments, we can assume without loss of generality that 
the fi converge to /o pointwise a.e. By continuous invertibility of r h- > the Uj converge to 

Mo pointwise a.e. Moreover, by construction, 

rb rb 



sup 



sup 



fi{x) dx < oo, 



so we can invoke Vitali's theorem 
convergence of ug to mq. 



recall the superlinear growth (57 1 



to conclude strong 
□ 



Proof of Proposition \2S\ It suffices to show that every subsequence of (uAj )jeN contains a sub- 
subsequence which is relatively compact. In view of Proposition |21[ we may thus assume — 
without loss of generality — that (wAj )jeN converges uniformly w.r.t. t € [0,T] in W2 to a curve 
€ C^/^([0, T]; W2). The verification of (65) then becomes an easy exercise, which is left to the 
reader. 



We verify ( 64 ) . As remarked in ( 66 ) , we have 



for all rt = 1, 



i?("A,)=Var^ {Hulf) 
Thus, the regularity estimate (l58|) implies 



with C given in (|59[), and with E, a from (61). This yields the uniform bound (64) 



Thus Theorern]23| applies and provides relative compactness of (mAjO^gn in L^{[Q,T] x /). 
Since L^-convergence implies weak convergence, it actually follows that ua^ converges to in 
L^([0,r] X /). Without loss of generality, we may even assume that u^. converges to a.e. on 
[0,r] X /. By continuity of P, also P('itAj) converges to P(m*) a.e. on [0,T] x /. Further, 



n=l 

< 2(6-a)C(E,a,T- 



M 

b — a 



Var^ (Hulf) 



which is j-uniformly bounded because of the regularity estimate ( |58| . By the growth property (57 1 
of <i>, we can invoke Vitali's theorem to conclude that P('Ua ) tends to P(m*) in L^{[Q,T] x /). □ 



22 



6.3. Weak formulation. Combining the compactness results from Proposition [2T] and Proposi- 
we know that every subsequence of (Aj)jgN contains a sub-subsequence for which UAj 

e C^/^[0,T];W2) n L\[0,T] X I), 



tion 



converges to some limit 
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uniformly w.r.t. t G [0,r] in W2, and strongly in L^{[0,T] x /); finally, also P(uAj) converges to 
P(m*) in L^([0,T] x /). To simplify notations, we denote that sub-subsequence simply by (ua), 
bearing in mind that A ~ Aj. In this subsection, we prove that every such limit is a weak 
solution to the initial value problem ([T]). 

Proposition 25. satisfies the weak formulation 

u^dt(pdxdt= / / u^Vxdx(p dx dt — / / P{u^)dxx'P dx dt (67) 



for all test functions ip from 

V:^ {ipe C°°([0,r] X /) I suppv? C (0,T) X /, ipx{t;a) = ipxit;b) = fa. t G [0,r]}. 
Moreover, attains the initial datum weakly-k as t 

Remark 26. Since weak solutions to ([T]) are unique (this follows, e.g., by metric contraction of 
the gradient flow), we conclude a posteriori that the entire sequence (uAj )jeN converges to the 
solution . 



Remark 27. By our definition (68 1 of test functions, the weak formulation (67) automatically 
induces homogeneous Neumann boundary conditions on u. 



The weak formulation ( 67 1 is obtained in the limit j — > 00 from a certain fully discrete variant 
of the weak formulation. The latter is derived by studying suitable variations of the minimizers 
u\_. In order to discuss this perturbation, let p G C°°{I) with Px{a) = Px{b) = be given, and 
let K > be such that 

\px{x)\<K, \pxx{x)\<Hi, \pxxx{x)\<K. for aU x G /. (69) 

Further, fix j G N and also n G N. We shall omit the index j in the following. 
Introduce the functionals A : X ^ M and : y ^ M by 

pM 

A(X)=/ p(X)dC and A^(x) = A(X5[x]). 
The derivative of A^ is given by 

pM 

[d^A^{l)\^= / p4X^^.[x])0,.dC for 

The variations we study are those induced by the gradient vector field v := W^^SijA^, i.e., 

Ww(x) = ajAg(x). (70) 
This vector field has a nice asymptotic expansion. 
Lemma 28. For every x G y, define pi(x) G M^^^ by 

Px{x) = {px{xi),px{x2), ■ . .,Px{xK-l))- (71) 
Then the residual vector v ~ w(x) — aI;(x) satisfies 

v'^Wv < CK^a{$,)MS{xf, (72) 



r-M 

k = l,...,K -I. 



with K defined in (69) and some universal constant C. 
Proof. We start by estimating 

p,:^wiy^ w(«(x) - p;[x]) =. axA^(x) - wp;[x]. 

By definition, we have 

X^[x](0 ^x, + ^^~/^-^ (e - Cfe) = + ^^^^(e - a-i) 
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for all ^ e [£,k-i,S.k\- Using the explicit form of W given in (27), we calculate: 



IJ'k 



2 1 
Px(X^[x]) - {-Px{xk) + -px{xk-i)) 



Okd^ 



r 2 1 1 

Pa;(X|[x]) - {-Px{xk) + 7^Px{xk+l))\dk 



{Xk - Xk-l) 5k 



1 ^ 



^Pxx (■^) 



Pxx{x) 



(1 - s)ds 



3 ' ' 3 

''^Pxx^^ ^ ^) ~h ^^^^^^ g 



(1 - s)ds, 



where x,x € I denote suitable intermediate values, depending on s. It follows that there is a 
universal constant Ci such that 

\pk\ < Cisup\pxx{x)\ max (4+4+i)^(x)^ 

m=l,...,K 



for every k — 1, . . . ,K — 1. Recalling the lower estimate on W in (28), it follows for z/ = W ^/i 
that 



proving our claim (|72|. 



6 



K-l 



2 ^ liCW 



k=l 



(4+1 + Skf)5ixf < 24C^,^ai^)MS{xf, 



Lemma 29. With k satisfying (69), the following estimate holds: 
^ (A^(xl) - A^(x'-i)) < -95cEc(51)V.(^l) 



+ ^(x2-x^-irw(x'A-x'A-i) 



□ 

(73) 
(74) 
(75) 



Proof. A Taylor expansion of p yields for arbitrary X, X' e X: 



A(X') > A(X) + / px{X) . (X' - X) dC - - / (X' - X)2 de 
JO ^ Jo 

With X = X|[x'A] and X' = X5[x;^"^], we obtain 

-(A«(xA)-A«(x'A-^)) -[xl-xr']„ / p.(X^[icA])0™d^ 

^ m=l -'O 



+ ^(^-^"'fw(xA-xr') 



Using the definition (70) of v, the symmetry W"^ = W, and the Euler-Lagrange equations (34), 

1 



the sum can be rewritten as 

A'-l 



(xl - xl-i)^Wz?(xl) = dsE^{x^^r'v{5r^). 



A) 



With the notations introduced in Lemma 28 we obtain further 

< dsE^{x'^fp-{x^^) + {dsE^{x^^fW-'dsE^{xl)f^{{,ylfW,,lf\ 



where the Cauchy-Schwarz inequality has been applied in the last step. The claim (73) now follows 
directly from the estimate (72). □ 
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Lemma 30. For every x e y^. 



P{u^[x]{x))pxxix)dx + / Vj:{x)px{x)u^[5(\{x)dx 



(76) 
(77) 



where € / are x-dependent quantities satisfying 

\x — x\, \x — x\ < (5(x). 



Proof. Relation ( 76 1 is a direct consequence of 



p{xk8k + Xk-idk-i) 



p{x) — 

Xk-i "^k Xk—1 



dx. 



which follows by a change of variables x = XkOk{C) + Xk-iOk-i{0- To prove (771, first observe 
that (30) implies 



K-l 



55fE^(x)^p;(x) - - ^ 



k=l 



■ Xk 



V 4+1 



Xk ~ Xk-l 



Px{xk) 



M 



K~l 



l^,(Xg[x]) ^p,(a;fe)0fede 



k = l 



We consider both terms on the right hand side separately. For the first, we obtain, using that 
-0'(f/r) = -P{r) and that Pxixi^) = Px{xk), 



K-l 



k=l L '^^^ 



Xk - Xk-1 

Sk 



K 

,(:r,)=^^'( 



Xk - Xk-l 



k = l 



_^Jxk-i ^Xk Xk — l 



{Px{xk) - Px{xk-l)) 



Pxx{xk)dx, 



with a suitable Xk G {xk-i, Xk) by the intermediate value theorem. For the other term, we perform 
a change of variables: 



Vx{xkQk + Xk-i0k-i){px{xk)9k+ Px{xk-i)Ok-i)d^= / Vx{x)px{x) 



dx. 



Xk ~ Xk-l 



with some x-dependent intermediate value x € {xk-i, Xk)- Summation over k — 1^ . . . ^ K provides 
([77|. □ 

Lemma 31. Let d g C^{0,T) be a non-negative test function of compact support in (0,T). Then 



J I 



■d'{t)p{x)u4t;x)dxdt < / i}{t)[- P{u^)pxx + VxPxu] dx dt. 

Ji 



(78) 



Proof. Multiply inequality ( 73 )-( 75 1 by r??(nr) > 0, and sum over n = l,...,Nr. On the left-hand 
A^(x^)-A^(x^-i) _ 



side, it follows by means of (76 1 that 



n=l 



N^-1 

•E 

n=0 
T 



■diin + 1)t) - ■&{nT) 



J I 



'd'^{t)p{x)uA{t; x) dx di. 



where the sequence of piecewise constant functions i?!^ converge to d' uniformly on [0,T]. The 
strong convergence of ua to u^, in is sufficient to pass to the limit j — >■ oo. 
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On the right-hand side, the first term can be rewritten using (77 1: 



Jo Ji 



where, for (n — l)r < t < nr, we know that |x — a;|, |a; — i;| < S{x^). The convergence (63 1 imphes 
that Pxxix) and px{x) converge to their respective Hmits pxxix) and Px{x) uniformly in a; e /. 
Likewise, the piecewise constant interpolants 7?r converge to d uniformly on [0,T]. The strong 
convergence of ma and of P{ua) to their respective limits and P(u*) in L^{[0,T] x /) thus 
suffices to pass to the limit with the integral. 

We still need to estimate the remainder terms, resulting from (74)&(75). Observe that 



A 



W(x 



A 



71=1 



1 ^ _ 

^^EW2(u«K],u«[x'-i])'<r(E-E) 



71=1 



by (45 1 and condition (61). Recalling that r depends on j, with tj l 0, this expression vanishes 
in the limit. Similarly, by (46) and (61), 



71=1 ^ n=l 

< (2r)i/2(E-E)i/^ 



1/2 



which is j-uniformly bounded. Hence, the remainder term resulting from (75 1 vanishes in the limit 
j — )• oo because of (63) and of (61). 



□ 



Proof of Proposition's^ First, observe that inequality (78) is actually an equality, since p can be 
replaced by —p everywhere. Next, recall that functions f of the type ip(t,x) = •d{t)p{x), where 
■d £ C^(0,r) is non- negative, and p € C°°{I) satisfies Px{a) = Px{b) = 0, are dense in the set V 
of test functions. Thus, the weak formulation (67) holds for all £ P. 

It remains to prove that the solution attains the initial datum vP weakly-* as t 0. However, 
this is a trivial consequence of the Holder regularity of € C^/^([0, T]; W2), of the uniform 
convergence of (uAj )jeN to u^, w.r.t. t £ [0,T] in W2, and of the approximation of vP by ■ I— ' 



7. Numerical results and proof of consistency 
7.1. Implementation. 

7.1.1. Choice of the initial condition. The numerical scheme is phrased in Lagrangian coordinates: 
the discretization £ = (^0, ^ii • • • 7 ^k) of the reference domain [0, M] is fixed, whereas the corre- 
sponding grid points x" — (a;", . . . ,x]^_i) S y on the interval / evolve in (discrete) time. In the 
numerical experiments that follows, our choice for the discretization of the initial condition is to 
use an equidistant grid x" with K vertices on /, 

xl ^ a + k{b - a)/K, 

and an accordingly adapted mesh ^ on [0, M], with 

^^^U°{xl), where U°(x)^( u°{y)dy for all a: € / 

J a 

is the initial datum's distribution function. This discretization has the property that 

u^{x)dx— / u%{x)dx for all fc = 1, . . . , if. 
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7.1.2. Time stepping. Each (time) step in the numerical scheme consists of solving the system 



(34 1 of Euler-Lagrange equations. In practice, this is done with a damped Newton method, which 



guarantees that the constraint ic^ € y — i.e., that a < a;" < • • • < x^j^_-^ < b — is propagated 
from the n — 1st to the nth iterate. To be more precise, recall that 



1 



d^E^.rix) = -W(x- jf^-^) + d^E^ix) 
is the functional whose unique root in y defines the nth time iterate x^, and that 

5|Ef ^(x) = -W + a|E.(x) 

T 

is its Jacobian. Given x^~^, we calculate x^ by means of the following algorithm: 

repeat 

dx:=-(a|E^.,(x))"'93cEg,,(x); 

while x+dx^y dx:=0.5dx; end; 

X := X + dx; 
until II dx||;i < tol and ||(35jE^.t-(x)||;i < tol; 
xl X. 

In our experiments, we use tol = 10^®. For the evaluation of 95jE^ ^(x) above, an explicit expression 
for the integrals 



voiL^ [x] dn = K o [x]e^ (0 dc 



is needed, see (30 1. Denoting by SO an anti-derivative of V , one finds 



o X^[x]0„(O d^ - / V^{x) dx 



XjYi ^rn— 1 



ViXrn) 



V{x^) - ^{x 

m— 1 J 
•^7n ^m — 1 



and analogously for the integral from to ^m+i- In combination, we obtain 



d,J Fox^[x](e)d^ 

Jo 



{Xk - .Tfe_l)2 



(2J(xfe)-53(a;fc_i)) + 



4-1 



{Xk+l - Xk} 



:{m{xk+i)-^{xk)) 



V{xk) 



Sk 



4+1 



^Xk — Xk-l Xk+l—Xf;^ 

A similar expression is available for the respective contribution to the Hessian i9|E^^^: 



-2V(xk) 



Sk 



Sk + i 



. ^ ( Sk 

(Xk-Xu-lY ' (Xk + l-Xk)^ J ' ^ y{Xk-Xk-l)^ (Xk + l-Xky J ^ 



S+ik 



m = k 



0, otherwise 
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Space X 




Space X 



Figure 2. Left: evolution of the (reference) solution ma with initial condition 
( 79| at times t = 0, r, . . . , 20r = 0.2, with time step r = lO'^ and K = 5000 grid 
points. The dotted line shows the stationary solution. Right: associated particle 
trajectories. 



7.2. Numerical experiments and convergence. The following numerical experiments are per- 
formed for the porous medium equation with quadratic nonlinearity, 

dtu = {u'^)xx + {VxU)x, 

on the interval / = [—1,1]. For the potential V, we choose 

V{x) — cos(7ra;), 

TT 

and as initial datum, we take the following function of unit mass M = 1: 

= C(-cos(27ra;) + 1.5) ((x + 0.5)4 + 1) with g ^ ^40 - 2807rJ^+ 423^^ ^ ^^^^ 

7.2.1. Reference Solution. Our numerical reference resolution is calculated with K — 5000 spatial 
grid points and a time step size r — 10~^. Figure [2]/left shows snapshots of the reference solution's 
spatial density after the first couple of time steps. One observes the typical behaviour for nonlinear 
drift diffusion equations: on a very short time scale, diffusion reduces the extrema of the initial mass 
distribution; subsequently, the drift dominates and transports the mass towards the equilibrium 
(dotted line) on a longer time scale. Figure [2]/ right displays the corresponding particle trajectories 
in the Lagrangian picture, i.e., how the points x]J move with (discrete) time n for fixed k. 

7.2.2. Fixed t. In a first series of experiments, wc fix the time step r — 10^^ and vary the 
number of spatial grid points K. In Figure [sj/left, the corresponding L^-distances to the reference 
solution ttref obtained in |7.2.1| above are shown as a function of time. Note that the counter- 
intuitive dramatic decay of the error for small times is explained by the strong contractivity of 
the nonlinear diffusion in in a neighborhood of the initial condition. Figure [s]/ right shows the 
L^-errors at T = 0.2. The observed convergence rate is of order K~^. 

7.2.3. Fixed parabolic mesh ratio. Next, we study the decay of the L^-error under mutual refine- 
ment of space and time. As it is standard in numerical experiments on parabolic equations, we 
fix the parabolic mesh ratio K^t. The value of this ratio is chosen such that the inverse CFL 
condition is satisfied in every experiment. In Figure [4j the error is plotted — similar as in the 
previous experiment — as a function of time (left) and at the fixed terminal time T ~ 0.2 (right), 
both for various choices of r. The observed order of convergence is y^r, which is in agreement 
with the result of experiment |7.2.2[ 
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Figure 3. Numerical error analysis with fixed time step r = 10 ^, using K = 
25, 50, 100, 200, 400, 800, 1600 grid points. Left: evolution of the i^-error \\uA{t)- 
Urof(i)|lLi(/)- Right: order of convergence at terminal time T = 0.2. 




Figure 4. Numerical error analysis with fixed parabolic mesh ratio K^t k, 0.257, 
using r = 5 • 10-^ 10-^ 5 • lO"", lO'^, 5 • lO'^, lO'^, 5 • 10"^, 10"!. Left: evolution 
of the i^-error ||ua(0 ^ Wrcf(i)|lLi(/)- Right: order of convergence at terminal 
time T = 0.2. 



7.2.4. Weakly convergent initial datum. In order to illustrate that it sufficies to approximate the 
original initial condition by its discretizations vP^ just weakly in L^{I), we use perturbed discrete 
initial data ^ that are biased by high-frequency oscillations of fixed amplitude 0.1, as indicated 
in Figure [5]/left. As expected, the perturbation becomes almost invisible already after the first time 
step, and the discrete solution ua,£ is indistinguishable from the one computed with unperturbed 
initial conditions u^. 

7.2.5. A discontinuous initial datum. For the last two series of experiments, we change the initial 
condition This first series is carried out with the discontinuous inital datum 

= if N > 0-75 or N < 0.25, ^^^^ 

0.9, otherwise. 



Similar to experiment 7.2.2 we fix r = 10 ^ and vary the number of grid points K. Figure [5]/ right 



displays the corresponding L^-error over the time interval t G [0,0.8]. In contrast to experiment 



7.2.2 the approximation error is zero initially, since the step function can be discretized exactly. 



However, the error jumps to a positive value (that is of the same order as the initial error in 
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Figure 6. The merely non-negative initial condition u*^ from (81) is approx- 
imated by strictly positive data -I- e. Left: discrete initial profiles for 
e — 10~^, 10"^, 10"'^, 10~^, 10~^. Right: qualitative behaviour of corresponding 
discrete solutions at T = 0.6, using r — 10^^, K = 200. 



experiment 7.2.21 in the first time step. Afterwards, the qualitative behaviour is very similar to 



that in experiment 7.2.2 The observed order of convergence (at T = 0.2) is again K ^. 



7.2.6. A non-positive initial datum. For this last series of experiment, we consider the initial 
condition 



u^\x) = ( - cos(27ra;) + 1.5) {{x + 0.5)^ + l) x 



-ix-0.5){x + 0.5) |x|<0.5 
|a;|>0.5 



(81) 



which vanishes outside of the subinterval [—0.5, 0.5] C /. The numerical scheme is not directly 
applicable to u^, but to any of its strictly positive approximations u° -|- e, see Figure p /left. The 
qualitative numerical results at T = 0.6 for various choices of e > are given in Figure 6]/right. 

7.3. Order of consistency. The experimental observations that the discrete solutions seem to 
approximate the reference solution with an error of order Kt can be supported theoretically by 
the following consistency consideration. 
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Assume that X : [0,T] x [0, A/] ^ / is a smooth solution of (10). Consider a discretization 
A = (t; ^) that is equidistant w.r.t. ^, i.e., = Mk/K with some ii' e N for all /c = 0, 1, 2, . . . , 
and consequently 5 := 5i = ■ ■ ■ = 5k- 

From X, we define a discrete "pseudo-solution" xa by restriction, i.e., :— X{nT,k5). We 
are going to show that xa satisfies the discrete evolution equation ( 34 1 up to an error of order 
6{0{t)+0{6'^)). To this end, we perform a Taylor expansion of X around a fixed point p :— {nT,kS) 
w.r.t. ^: 

x^i = X(p) ± 5X^ip) + yXee(p) ± jX^i^ip) + 0(5% 
For the diffusion term, we find 



[^r{Mp))Xmip) + l^"'iMp)Wp)') + 



52 /I 



so that 



For the drift term, we obtain 



f^'^'^' y4X^[xA](e))efc(e)de - ±S C V,{{1 - s)xl + - ds 



±5 / [F,(X(p)) ± sJF,,(X(p))X5(p) + 0(5^)] (1 _ s) ds 
Jo 

±^F,(X(p)) + ^F,,(X(p))xap) + 0{5^), 



so that 

f.(fc+l)<5 



In combination, 



Moreover, we have 



and likewise 



y4Xg[xA](0)^fc(0 = 5v,{np)) + (^(<5'). 

(fe-l)(5 



[a^jE^Cx)]^ -<5(^'(X5(p))^ + F,(X(p)) + 0(<52)). (82) 



[W^'^\ = 5X{p')+0{5^), 
where p' — {{n — 1)t, kS). Finally, using 

X{p)~X{p')=TXt{p)+0{T'), 

we obtain the relation 

1 



- [W(x'A - ^^-')] , = S{X,{p) + 0{5^) + 0(r)). (83) 



r 



Using the continuous evolution equation ( |10| ) in ([83| and ( [82| leads to 

^ [W(x'A - ^^-i)] , - - [d^M^ + '5(0(r) + 0{5^)). 
for all admissible /c and n. 
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